\chapter{Carroll2006}
\section{问题的描述}
经济学中最优化问题通常表述为，$ \max\quad \sum_{s=t}^{T} \beta^{s-t}u(C_s)$,如果把效用函数用常风险回避形式，则可以写为，
\[ \max\quad \sum_{s=t}^{T} \beta^{s-t}\frac{C^{1-\rho}}{1-\rho} \]

这样的表达在微观中，往往是异质性的劳动力冲击，宏观中则是代表性行为人的加总生产率冲击。\textbf{期末}的资产是市场资源减去消费，即
\begin{equation}\label{eq_a}
A_t=M_t-C_t
\end{equation}

这个资产在下一期\textbf{开始}的时候会有折旧，在宏观中有折旧，在微观中通常没有折旧：
\[ K_{t+1}=A_t(1-\rho) \]

然后从 $ t+1 $期开始，此时$ K_{t+1} $还未用来生产产出，到这一期中间的时候，产出生产出来，然后就可以合并到市场资源中，即，
\begin{equation}\label{eq_m}
M_{t+1}=\underbrace{\Theta_{t+1}P_{t+1}}_{\equiv L_{t+1}}W_{t+1}+K_{t+1}R_{t+1}	
\end{equation}

这里，$\Theta_{t+1} $是i.i.d.的暂时性冲击(比如失业)，已被标准化满足$ E_t\Theta_{t+n}=1,\forall n>0 $，而在宏观中，它就等于1。永久性劳动生产率$ P_{t+1} $如下演进，
\[ P_{t+1}=\Gamma_{t+1}P_t\Psi_{t+1} \]
永久性冲击$ E_t\Psi_{t+n}=1\forall n>0 $，$ \Gamma_t $外生可预期。

通过除以永久性劳动生产率可以将\eqref{eq_a}式和\eqref{eq_m}式标准化，即
\begin{align}
	a_t=&m_t-c_t\\\label{eq_mtp1}
	m_{t+1}=&\underbrace{\Theta_{t+1}}_{\equiv l_{t+1}}W_{t+1}+\underbrace{\frac{a_t(1-\rho)}{\Gamma_{t+1}\Psi_{t+1}}}_{k_{t+1}}R_{t+1}
\end{align}

通过这些表述，可以把最优化问题归纳为一个贝尔曼方程，
\[ V_t(M_t,P_t)=\max_{C_t}\{u(C_t)+\beta E_t[V_{t+1}(M_{t+1},P_{t+1})]\} \]
或者写成如下形式\footnote{注意值函数是$ 1-\rho $次齐次的。}，
\begin{equation}\label{eq_v}
v_t(m_t)=\max_{c_t}\left\{u(c_t)+\beta E_t\left[(\Gamma_{t+1}\Psi_{t+1})^{1-\rho}\cdot v_{t+1}\left(\underbrace{l_{t+1}W_{t+1}+R_{t+1}\frac{a_t(1-\rho)}{\Gamma_{t+1}\Psi_{t+1}}}_{=m_{t+1}}\right)\right] \right\}	
\end{equation}

现在可以定义一个新的变量——\textbf{后一期值函数}，为，
\[ \mathcal{V}_t(a_t)=\beta E_t\left[(\Gamma_{t+1}\Psi_{t+1})^{1-\rho}\cdot v_{t+1}\left(l_{t+1}W_{t+1}+R_{t+1}\frac{a_t(1-\rho)}{\Gamma_{t+1}\Psi_{t+1}}\right)\right]  \]
于是\eqref{eq_v}式可以写为，
\[ v_t(m_t)=\max_{a_t}[u(m_t-a_t)+\mathcal{V}_t(a_t)] \]
注意到$ \mathcal{V}_t $对$ a $的偏导$ \mathcal{V}^a_t $为，
\begin{align*}
	\mathcal{V}^a_t=&\beta E_t\left[(\Gamma_{t+1}\Psi_{t+1})^{1-\rho}\cdot v_{t+1}^m(m_{t+1})\cdot \frac{R_{t+1}(1-\rho)}{\Gamma_{t+1}\Psi_{t+1}}\right]\\
	=&(1-\rho)\beta E_t[(\Gamma_{t+1}\Psi_{t+1})^{-\rho}v_{t+1}^m\cdot R_{t+1}]
\end{align*}
包络定理意味着,
\[ v_t^m(m_t)=u'(c_t) \]
那么一阶条件欧拉方程为，
\begin{align}
u'(c_t)=&\mathcal{V}^a_t(a_t)\\
=&(1-\rho)\beta E_t[(\Gamma_{t+1}\Psi_{t+1})^{-\rho}u'(c_{t+1})\cdot R_{t+1}]	\\\label{eq_euler}
=&(1-\rho)\beta E_t[u'(c_{t+1}\Gamma_{t+1}\Psi_{t+1})\cdot R_{t+1}]\qquad\text{因为消费函数的设置}
\end{align}
\section{解法}
\textbf{标准的解法}，就是假设$ c_{t+1} $已知，然后设置一些有序格点$ \mu_i\in \vec{\mu}\equiv \{\mu_1,\mu_2,\cdots,\mu_I\} $，这些格点代表的是状态变量$ m_t $。然后，对于每一个格点，对\eqref{eq_euler}式数值求根，
\begin{equation}\label{eq_chi}
	u'(\chi_i)=\mathcal{V}_t^a(\mu_i-\chi_i)
\end{equation}
点对$ \{\mu_i,\chi_i\} $就可以用来构造插值以逼近$ c_t $。问题在于这种解法在每个特定格点上都要解\eqref{eq_euler}式。
\begin{table}[H]\centering
\begin{tabular}{cc}\hline
变量&格点\\\hline
$ m_t $& $ \mu_i $\\
$ c_t $&$ \chi_i $\\\hline
\end{tabular}
\end{table}

\textbf{内生格点法}可以避免每次迭代要对每个点求解\eqref{eq_euler}式。他的方法在于不是直接格点化$ m_t $，而是格点化$ a_t $为$ \alpha_i\in \vec{\alpha}\equiv\{\alpha_1,\alpha_2,\cdots,\alpha_I\} $，而对于每个格点的\textbf{后一期值函数}的导数$ \mathcal{V}_t^a(\alpha) $是很容易计算的，算出来后，对它求逆就能得到$ \chi_i $，
\[ \chi_i =u'^{-1}[\mathcal{V}_t^a(\alpha_i)] \]
利用下式就能得到$ \mu_i $，
\[ \mu_i=\alpha_i+\chi_i \]
于是我们再次拥有了$ \{\mu_i,\chi_i\} $对，从而又可以插值。具体算法我们总结如下，
\begin{enumerate}
	\item 格点化$ a_t $，然后初始化$ m_t,c_t $，记为$ m_t^0,c_t^0 $，比如$ m_t^0,c_t^0 $分别为一个包含20个元素的向量，值均为0，20即为对应$ a_t $的格点数目。
	\item 使用\eqref{eq_euler}式的右边计算$ \mathcal{V}_t^a $。而\eqref{eq_euler}式的计算涉及到$ c_{t+1} $，此时利用初始化的$ \{m_t^0,c_t^0\} $对，注意到\eqref{eq_mtp1}式可以计算$ m_{t+1} $，然后插值获得$ c_{t+1} $，这样根据\eqref{eq_euler}式就能求得新的$ c_t $，即$ c_t^1 $，然后将此$ c_t^1 $替换掉$ c_t^0 $。
	\item 利用$ m_t^1=a_t+c_t^1 $，更新$ m_t^0 $，至此获得新的$ \{m_t^1,c_t^1\} $对。
	\item 重复第2和第3步，直至收敛。
\end{enumerate}

\section{R代码}
以下为我对Carroll(2006)matlab代码修改后的R代码：
\begin{lstlisting}[language=R,numbers=left,firstnumber=33,mathescape=false]
rm(list = ls())
library(reshape2)
library(ggplot2)
library(magrittr)
library(stringr)
devtools::load_all()

# ---------函数--------
# GothVP is defined as 后一期值函数的导数
GothVP <- function(a, tranval,permval,Dep,Beta,R,W,M,cons){
	
	# EUP 是效用函数导数的期望
	EUP <- matlab::zeros(matlab::size(a))
	
	# 这里的循环主要是为了计算期望，需要遍历各个概率
	for (i in 1:length(tranval)){
		for (j in 1:length(permval)){
			CapPsi <- permval[j]
			CapXi <- tranval[i]
			k <- a*Dep/(G*CapPsi)
			EUP <- EUP+Dep*Beta*R(k,CurlyEpsilon)*
			UP(G*CapPsi* signal::interp1(M,cons,R(k,CurlyEpsilon)*k + W(k,CurlyEpsilon)*CapXi,extrap = T),Rho)*
			tranprob[i]*permprob[j]
		}
	}
	return(EUP)
}

NP <- function(cc, Rho) cc^(-1/Rho) # u'的反函数
UP <- function(cc,Rho) cc^(-Rho) # u'

# --------求解----------
# 参数设置
Beta <- 0.96   # time preference
Rho <- 2       # coefficient of relative risk aversion
n <- 20        # number of grid points
G<-1.01      # permanent income growth rate
Dep<-0.90    # depreciation rate
CurlyEpsilon<-0.36  # Cobb-Douglas production function parameter
permval<-c(0.90,1.00,1.10)   # permanent shock values
permprob<-c(.25,0.50,0.25)  # permanent shock probabilities
tranval<-1.00             # transitory  shock values
tranprob<-1.00            # transitoty shock probabilities
R <-function(k,CurlyEpsilon) 1+CurlyEpsilon*k^(CurlyEpsilon-1)   # set up interest rate function
W <- function(k,CurlyEpsilon) (1-CurlyEpsilon)*(k^CurlyEpsilon)   # set up wage rate function
kSS <- ((((G^Rho)/(Beta*Dep))-1)/CurlyEpsilon)^(1/(CurlyEpsilon-1))
aSS <- kSS*G/Dep
AlphaVec <- exp(seq(0,log(3*aSS),length.out = 20))-1       # set up grid points
PeriodsToAdd=99 # 迭代期数

# 初始化
cons <- seq(0,n-1) %>% matrix(ncol = 1)
M <- seq(0,n-1) %>% matrix(ncol = 1)

# 求解
for (i in 1:PeriodsToAdd){
	# calculate ct from each grid point in AlphaVec
	ChiVec <- NP(GothVP(AlphaVec,tranval,permval,Dep,Beta,R,W,M[,ncol(M)],cons[,ncol(cons)]),Rho) # inverse Euler equation
	MuVec  <- AlphaVec+ChiVec
	M <- cbind(M, matrix(MuVec,ncol = 1))             # matrix of interpolation data
	cons <- cbind(cons, matrix(ChiVec,ncol = 1))             # matrix of interpolation data
}
\end{lstlisting}

\section{另一个简单例子}
一个RC经济，通常的约束条件为，
\begin{equation}\label{eq_ma}
\underbrace{k_{t+1}}_{a_t}=\underbrace{f(k_t)+(1-\delta)k_t}_{m_t}-c_t	
\end{equation}


一般的贝尔曼方程写为，
\[ V(k_t)=\max_{k_{t+1}}[u(c_t)+\beta V(k_{t+1})] \]
也可以写成，
\[ V(m_t)=\max_{a_{t}}[u(m_t-a_t)+\beta V(m_{t+1})] \]

%然后定义\textbf{后一期值函数}为，
%\[ \mathcal{V}(a_{t+1})= \beta V(m_{t+1})\]
一阶条件为，
\begin{align}
u'(c_t)=&\beta V_{m_{t+1}}(m_{t+1})\cdot\frac{\partial m_{t+1}}{\partial a_{t}}\\\label{eq_key}
=&	\beta V_{m_{t+1}}(m_{t+1})\cdot [f'(k_{t+1})+(1-\delta)]
\end{align}

而包络定理表明，
\[ V_{m_t}=u'(c_t) \]
因此，\eqref{eq_key}式可以最终写为(也就是欧拉方程)，
\begin{equation}\label{eq_el2}
u'(c_t) =\beta u'(c_{t+1})\cdot [f'(k_{t+1})+(1-\delta)]	
\end{equation}

\textbf{标准的解法}就是，
\begin{enumerate}
	\item 格点化$ k_t $，猜测一个$ c_t $序列，形成数值对序列$ \{k_t,c_t^0\} $
	\item 基于\eqref{eq_ma}式，利用$ c_t^0 $,计算得到$ k_{t+1}^0 $
	\item 利用数值对$ \{k_t,c_t^0\} $和$ k_{t+1}^0 $插值得到$ c_{t+1}^0 $
	\item 于是解\eqref{eq_el2}式得到$ c_t^1 $，更新$ c_t^0 $
	\item 重复第2步至第4步，直至收敛
\end{enumerate}

\textbf{内生格点法}就是，
\begin{enumerate}
	\item 格点化$ a_t $，猜测一个数值对序列$ \{m_t^0,c_t^0\} $
	\item 因为$ a_t $就是$ k_{t+1} $，因此，利用$ a_t $可以得到$ m_{t+1}^{0} $，然后基于数值对$ \{m_t^0,c_t^0\} $插值得到$ c_{t+1}^0 $
	\item 利用$ a_t,c_{t+1}^0 $和\eqref{eq_el2}式得到$ c_t^1 $，从而更新了$ c_t^0 $
	\item 利用$ m_t^1=a_t+c_t^1 $更新$ m_t^0 $
	\item 重复第2步至第4步，直至收敛
\end{enumerate}

\textbf{可以看到两种方法真正的区别在于因为$ a_t $的设置，使得期望的计算大为简化，要知道如果\eqref{eq_el2}式中涉及到t+1的脚标都是要计算期望的\footnote{我这个例子用的是确定性模型，随机模型是必须要计算期望的。}，标准解法必须计算所有格点上的期望，而内生格点法只需要在插值$ c_{t+1} $的时候计算一下。}

以后这样的最优规划，只要写出欧拉方程\eqref{eq_el2}式以及约束方程\eqref{eq_ma}式，对照此步骤即可。

基于RBCofABC第六章的RBC模型，我写了一个内生格点求解的R代码(\lstinline|EGM.R|)如下，
\begin{lstlisting}[language=R,numbers=left,firstnumber=33,mathescape=false]
rm(list = ls())
library(ggplot2)
library(magrittr)
library(signal)
devtools::load_all()

# 参数设置
a <- seq(0.06,6,0.06)
beta <- .98
delta <-  .1
theta <- .36
numits <- 100

# 函数定义: u = log(cc), 对数效用下偏导和偏导的逆是一样的
NP <- UP <- function(cc) 1/cc

# 初始化
m <- cc <- seq(1,length(a)) %>% matrix(ncol = 1)

for (i in 1:numits) {
	# 插值计算c(t+1),10式中mt的定义
	m_next <- a^theta + (1-delta)*a
	c_next <- interp1(m[,ncol(m)],cc[,ncol(cc)],m_next,extrap = T)
	
	# 笔记中的13式反算出ct
	cc <- NP(beta * (theta * a^(theta-1)+1-delta) * UP(c_next)) %>% cbind(cc,.)
	
	# 更新m
	m <- (a + cc[,ncol(cc)]) %>% cbind(m,.)
}

# 通过m反算出k
mfun <- function(kt,mt) kt^theta+(1-delta)*kt-mt
k <- numeric(nrow(m))
for (i in 1:100) {
	k[i] <- uniroot(mfun,c(0,10),mt = m[i,ncol(m)])$root
}

# 画图
ggplot(data.frame(kt = k, knext = a),aes(x = kt,y = knext)) + geom_line()
\end{lstlisting}